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Abstract. Numerical methods are developed to simulate the wave propagation in het- 
erogeneous 2D fluid / poroelastic media. Wave propagation is described by the usual 
^ I acoustics equations (in the fluid medium) and by the low-frequency Blot's equations 

(in the porous medium). Interface conditions are introduced to model various hy- 
draulic contacts between the two media: open pores, sealed pores, and imperfect 
pores. Well-possedness of the initial-boundary value problem is proven. Cartesian 
grid numerical methods previously developed in porous heterogeneous media are 
adapted to the present context: a fourth-order ADER scheme with Strang splitting 
for time-marching; a space-time mesh-refinement to capture the slow compressional 
^ ■ wave predicted by Blot's theory; and an immersed interface method to discretize the 

interface conditions and to introduce a subcell resolution. Numerical experiments and 
comparisons with exact solutions are proposed for the three types of interface condi- 
tions, demonstrating the accuracy of the approach. 
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1 Introduction 

The theory developed by Blot in 1956 (SHH is largely used to describe the wave propaga- 
tion in poroelastic media. Three kinds of waves are predicted: the usual shear wave and 
"fast" compressional wave (as in elastodynamics), and an additional "slow" compres- 
sional wave observed experimentally in 1981 Il33ll . This slow wave is a static mode below 
a critical frequency, depending on the viscosity of the saturating fluid. In the current 
study, we will focus on this low-frequency range. 
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The coupling between acoustic and poroelastic media is of high interest in many ap- 
plications: sea bottom in underwater acoustics ||38H , borehole logging in civil engineer- 
ing II36J , and bones in biomechanics f:21]. Many theoretical efforts have dealt with the 
acoustic / porous wave propagation. Various boundary conditions have been proposed 
to describe the hydraulic contacts: open pores, sealed pores, and imperfect pores involv- 
ing the hydraulic permeability of the interface |l6l|15H36l. Reflection and transmission 
coefficients of plane waves have been derived 1139 II . The influence of the interface condi- 
tions on the existence of surface waves has been investigated in the case of inviscid IITSlI 
and viscous saturating fluids IITSIITTII in the porous material. The time-domain Green's 
function has been computed by the Cagniard-de Hoop's method 111211161 . Experimental 
works have shown the crucial importance of hydraulic contact on the generation of slow 
compressional wave Il35ll . 

The literature dedicated to numerical methods for porous wave propagation is large: 
see 1123 II . Il9ll for a review, and the introduction of MTTH for a list of time-domain methods. 
Coupled fluid / porous configurations have been addressed by an integral method IITSU , 
a spectral-element method 1132 II , and a pseudospectral method [37], to cite a few. To sim- 
ulate efficiently wave propagation in fluid / porous media, numerical methods must 
overcome the following difficulties: 

• in the low-frequency range, the slow compressional wave is a diffusive-like solu- 
tion, and the evolution equations become stiff |l34ll . It drastically restricts the stabil- 
ity condition of any explicit method; 

• the diffusive slow compressional wave remains localized near the interfaces. Cap- 
turing this wave - that plays a key role on the balance equations - requires a very 
fine spatial mesh; 

• an accurate description of arbitrary-shaped geometries with various interface con- 
ditions is crucial. These properties are badly discretized by finite-difference meth- 
ods on Cartesian grids. Alternatively, unstructured meshes provide accurate de- 
scriptions, but the computational effort greatly increases. 

• an accurate modeling of the hydraulic contact at the interface is also required. In 
particular, as far as we know, imperfect pore conditions still have not been ad- 
dressed in numerical models. 

To overcome these difficulties, we adapt a methodology previously developed in porous 
/ porous media 111.1 and fluid / viscoelastic media li28J . Three Cartesian grid numerical 
methods are put together. A fourth-order ADER scheme with Strang splitting is used to 
integrate the evolution equations, ensuring an optimal CFL condition of stability. Spe- 
cific solvers are used in the fluid medium and in the porous medium. Their coupling is 
ensured by an immersed interface method, that discretizes the interface conditions and 
provides a subcell resolution of the geometries. Lastly, a space-time mesh refinement 
aroimd the interfaces captures the small scale of the slow waves. 
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The article is organized as follows. In section |2l acoustics and poroelastic equa- 
tions are recalled. We introduce the interface conditions, and we prove that the initial 
boundary-value problem is well-posed. In section |3l numerical tools are presented. In 
sectionHl numerical experiments are proposed, based on realistic sets of physical param- 
eters. Comparisons with analytical solutions demonstrate the accuracy of our approach. 
In section|5l future lines of investigation are proposed. 



2 Physical modeling 

2.1 Acoustics and Biot's equations 




n 



Figure 1: Interface F separating a fluid medium Qg and a poroelastic medium Qj. 

Let us consider a 2D domain with a fluid medium Qq and a poroelastic medium Qi 
(figure [T]). The interface T separating Qq and Qi is described by a parametric equation 
{x{T),y{T)) (figure [U. Tangential vector t and normal vector n are defined at each point 
P along r by: 

i={x,y)\ n=(y',-xy. (2.1) 

The derivatives x' = || and y' = ^ are assumed to be continuous everywhere along T, 
and to be differentiable as many times as required further. 

In the fluid domain Qq, the physical parameters are the density p f and the celerity of 
acoustic waves c. The acoustics equations write 

3v „ „ 

g (2.2) 

where v = (^1,^2)^ is the acoustic velocity and p the the acoustic pressure; fp represents 
an external source term. 

The poroelastic medium Oi is modeled by the low-frequency Biot equations ||6l where 
the physical parameters are 
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• the dynamic viscosity rj and the density Pf of the saturating fluid. The latter is 
assumed to be the same than in Qq, hence the notation is used in both cases; 

• the density ps and the shear modulus }i of the elastic skeleton; 

• the porosity <(p <1, the tortuosity a>l, the absolute permeability k, the Lame 
coefficient of the saturated matrix A f, and the two Biot's coefficients /S and m of the 
isotropic matrix. 

The conservation of momentum and the constitutive laws yield 

V.£r = 0, 



(2.3) 



9vs 


9w 










Pf dt 





where Vs=^ 



-w+V 

K 

cr = Ce{xXs)-^pI, 
p = -m(^V.Us + V.W), 
= (usi/'^s2)^ is the elastic velocity, w=(p{vf—Vs 



dt 



- {wi, W2Y is the filtra- 
tion velocity, is the fluid velocity, cr is the elastic stress tensor, £(us) = \ (Vus+^Vug) 
is the elastic strain tensor, and p is the pressure. The following notations have also been 
used in (|Z3l): p,u = ^Pf, p = (ppf + {l-(p)ps, and 



/ Ao+2^ Ao 
2^ 
V Ao Ao+2^ / 



(2.4) 



where Aq = A^ — ra is the Lame coefficient of the dry matrix. 

To be valid, the second equation of (12.311 requires that the spectrum of the waves lies 
mainly in the low-frequency range, involving frequencies lower than 



Inaxpf 

If f>fc, more sophisticated models are required 



(2.5) 



that are not addressed here. 



2.2 Evolution equations 

A velocity-stress formulation of the evolution equations is obtained by differentiating the 
last two equations in (|2.3|) in terms of time t. Setting 

{ivi,V2,pf inQo, 
(2-6) 
{Vsl,Vs2,ZOi,ZV2,Crn,0-i2,Cr22,p) in Oi, 
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where an, ai2, and are the independent components of the stress tensor tr, one de- 
duces from (I2.2I) and (I2.3|) the first-order linear system of partial differential equations 
with source term 

3 3 3 

^U+A — U+B — U = -SU+F. (2.7) 

dt dx dy 

The system (12 .ZD is completed by initial values and radiation conditions at infinity. 

In \2.7\ , A, B and S are 3x3 matrices in Qq, and 8x8 matrices in Qi; the vector F 
accounts for the acoustic source in (12.211 . In Oq, S = 0, while in Oi the spectral radius of S 




f (Hz) f (Hz) 



Figure 2: Phase velocities (a) and attenuations (b) of Biot's waves, using the parameters given in table[T] p/: 
fast compressional wave; ps: slow compressional wave; s: shear wave. In (a), the horizontal dotted lines refer 
to the eigenvalues c^, c^g and cf of A and B. 

The non-null eigenvalues of A and B are real: ±c (acoustic wave) in Oq; (fast 
compressional wave), ±c^ (slow compressional wave), and ±cf (shear wave) in Oi, 
satisfying 0< max (c~,c^) <c^f These eigenvalues in Qi are the high-frequency limits of 
the phase velocities of the poroelastic waves: 

lim Cp;(/)=c~^, lim Cp,(/)=c;, limCs(/) = C. (2-9) 

The fast compressional wave and the shear wave are almost non-dispersive and non- 
diffusive solutions. On the contrary, the phase velocity of the slow compressional wave 
tends to zero with frequency (figure |2]-a); at higher frequencies, this wave propagates, 
but it is highly attenuated (figure |2}-b). For a detailed dispersion analysis, the reader is 
referred to standard texts IHIHl. 
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2.3 Hydraulic contacts 

Four waves are involved in the acoustic / porous configuration: one acoustic wave in CIq, 
and three poroelastic waves in Qi . Consequently, four independent interface conditions 
need to be defined along T. Authors |l6Hl5l[l9H36l have proposed the following general 
conditions: 

' vo.n = Vsi.n+wi.n, (2.10a) 
— pon = (ri.n, (2.10b) 

[P] = -Y^' (2.10c) 
V |n| 

where the subscripts and 1 refer to the traces on the Qq or Oi sides, and [p] denotes the 
jump of p from Qq to Qi . The first scalar equation (I2.10a|) follows from the conservation 
of fluid mass. The second vectorial equation (|2.10b|l expresses the continuity of normal 
efforts. The third scalar equation (l2.10cD is a local Darcy's law. It models the hydraulic 
contact between the fluid and the porous medium, and involves an additional parameter 
/C, called the hydraulic permeability of the interface. The division by |n| ensures that 
the hydraulic contact is independent from the choice of the parametric equation of F. 
According to the value of /C, various limit-cases are encountered: 

• if /C ^ +0O, then the last equation in (|2.10|) becomes [p] = 0, modeling the commonly 
used open pores; 

• if /C ^ 0, then no fluid exchange occurs across F, and the last equation in (12.1011 is 
replaced by wi.n = 0, modeling sealed pores; 

• if < /C < +00, then an intermediate state between open pores and sealed pores is 
reached, modeling imperfect pores. 

The following proposition states that the interface conditions (12.101) coupled with the 
evolution equations (|2.2|) and (|2.3P yield a well-posed problem. 

Proposition 2.1. Let 

E = Ei + E2 + E3, 

with 




(2.11) 



7 



Then, E is an energy which satisfies 

'4 = -f !L„-in-fl.^^iY<o. (2.12) 

at Jci^K JtIC \n\ 

Proof. The first equation of (I2.2D is multiplied by v and integrated over Oq: 

The first term in (12.1311 writes 

By integration by part, and using the second equation of (|2.2|) , one obtains 
/ vVpdO = / vo-npoi^r— / V.vpcfQ, 

= /^vo.npodr+/^^-l,p|dn, (2.15) 
= /^vo.npodr+Al|^^_l_p2,n, 

which concludes the energy analysis in the fluid domain. The first equation of (|2.3P is 
multiplied by and is integrated over Qi: 

The first term in (|2.16ll writes 

By integration by part, and using the third equation of (|2.3|) . one obtains 

— / VsV.crdD. = / Vsi.(ri.ndT+ / e{vs):(rdO., 
Jcii Jr JQi 

= fvsi.(ri.ndT+ f e{ws):{CE{us)-^pl)dn, 

Jt Ja.1 

= [ Vsi.(ri.ndY+ [ ^e{us):Cs{us)dn- [ ^pV.VsdCl, 

Jt Jcii at Jcii 



jvsi.cri.ndl+^^ CE{us):E{us)dnj-J^ ^pV.v,-dn. 

(2.18) 
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The second equation of (|2.3|) is multiplied by w, and then it is integrated over Oi : 

L (p/w^+p.w^ + f w^+wVp) dn = 0. (2.19) 
The second term in (I2.19|) writes 

/ p^aw^dn = ^l f p^,vf^dn. (2.20) 

JDi at atlJcii 

By integration by part, and using the fourth equation of (|2.3|) . one obtains 

/ wVpdCl = — wi.npidF— / pV.wdO, 

Jcti Jt Jfli 

= -/^w,.np,dr + /^_4(ip + pV.u,)rfn, (2.21) 
When adding \1.\6) and (|2.19|) , it remains 
Equations (I2.13N2.22D provide 

^ = - / Iw^dn-it, (2.23) 
where xp is simplified from the conditions l|2.10|l : 

^ = ^(vo.npo+Vsi.ncrin-wi.npi)dr, 
= ^(vo.npo-Vsi.npo-wi.npi)dr, 

= ^((vsi.n+wi.n)po-Vsi.npo-wi.npi)(ir, (2.24) 
wi.n[p]dr. 



/C |n| 

which gives H2.12II . It remains to prove that E is a positive definite quadratic form. This 
is obviously the case of Ei and in E3. Definite positivity of C£(us) : £(us) is assumed 
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by hypothesis on the elastic tensor C IfMll . which concludes the analysis of £3. Lastly, the 
matrix 

9 Pf 
Pf Pw 



M 



is definite positive: 



detM = ppjv-pj, 



apfps>0, 



which concludes the analysis of £2- 



□ 



Each term in (12.111) has a clear physical significance: Ei is the acoustical energy, E2 
is the poroelastic kinetic energy, and E3 is the poroelastic potential energy; E3 is easily 
computed from (I2.4I) using the closed-form expression: 



C£(Us):£(Us) 



Ao 



■(C7ii+j6p)(c722+|6p). 



(2.25) 



2^{Ao+^) 

The decrease rate of the total energy is governed as usual by the intrinsic attenuation 
due to the viscous saturating fluid J^^ ^ w^dO, but also by the imperfect pore condition 

It X (^1 -n)^ (ir. In particular, even if the viscous effects are neglected {}] = 0), a part of the 
mechanical energy is dissipated by the interface in the case of imperfect pores. 

2.4 Matrix formulation of the interface conditions 

The various limit-cases in (|2.10|) are written in an alternative way, well-suited for further 
discretization of interface conditions (section 1331 and appendix |B]). First, the open pore 
conditions yield 

vo.n = Vsi.n+wi.n, 



-pon^ = (c^i-n).n, 

(£ri.n).t = 0, 

(£ri.n).n+pin^ = 0. 
Second, the sealed pore conditions yield 

vo.n = Vsi.n, 

-pon^ = (£^i-n).n, 

(£ri.n).t = 0, 

wi.n = 0. 



(2.26) 



(2.27) 
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Third and last, the imperfect pore conditions yield 



vo.n = Vsi.n+wi.n, 
-pon^ = (c^i.n).n, 
(£ri.n).t = 0, 

(tri. n). n + pi n^ + — wi.nlnl =0. 



(2.28) 



In the three cases (|2.26|) - !|2.28|l , the four scalar interface conditions are written as two jump 
conditions and two boundary conditions on the Oi side. 

These conditions are now written in a matrix way, useful for the further derivation of 
r-th order interface conditions (r > 1). On the side Q; (/ = 0, 1), the boundary values of the 
spatial derivatives of U up to the r-th order are put in a vector U-: 



lJ'i= lim 



-U^ 



(2.29) 



where Z = 0,...,r and m = 0,...,l. The vector U- has Wi, = 3(r+l) (r+2)/2 components in 
the fluid domain Qq, and rij, = 4(r+l) {r+2) components in the poroelastic domain Oi. 
Based on this formalism, the zero-th order interface conditions (I2.26D - II2.28D are written 



POttO_pUttU 



OttO 



L?U?: 



0, 



(2.30) 



where Cq is a 2 x 3 matrix, and Cj, L^j* are 2x8 matrices; they are detailed in appendix lAl 



3 Numerical modeling 

3.1 Integration of evolution equations 

The system (I2.7D is solved on a uniform Cartesian grid, with spatial mesh sizes Ax = Ay 
and a time step At. Due to the source term in the poroelastic medium Oi, a straightfor- 
ward discretization of (I2.7D is inefficient: a Von-Neumann analysis of stability gives 

(Ax 2 \ 

At<min ^, , (3.1) 

\max(cj -R(Sj / 



where the spectral radius R(S) can become large (12. 8D . It is more efficient to split (I2.7[l 
and to solve alternatively 

3 3 3 

^U + A — U+B — U = 0, (3.2a) 
dt dx dy 

~\ 

— U = -SU. (3.2b) 

at 
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The operators used to solve the propagative part (|3.2a|) and the diffusive part (|3.2b|l are 
denoted by and H^, respectively. Second-order Strang splitting is used Il24ll . hence 
time-marching can be written 

. u[;) = H,(f)u;;^, 

. = H«(A0U;]^ (3.3) 

. Uf;' = H,(f)ug^ 

The discrete operator Ha in (13.311 is a fourth-order ADER scheme Il22l . On a Cartesian 
grid, this explicit scheme amounts to a fourth-order Lax-Wendroff scheme. It satisfies 
the stability condition max(c) ^ < 1- The diffusive operator is exact and unconditionally 
stable: Hb(i) =e^^', see equation (18) in [11]. If dissipation processes are neglected, then 
Hfc is also the identity: it is always the case in Qq, and it is also true in Qi when ?/ = 0. 

When S 7^ (i.e. when t] ^ 0), the coupling between parts (I3.2all and (I3.2bll decreases 
formally the convergence rate from 4 to 2. In counterpart, the optimal condition of stabil- 
ity is recovered: max(c) ^ < 1/ which is essential for real parameters simulations. With- 
out splitting, the time step would be divided by about 3 and 4 in the Test 2 and Test 3 of 
section m respectively. 



3.2 Mesh refinement 

In the low-frequency range, the slow compressional wave behaves like a diffusive non- 
propagating solution, with very small wavelength (section |Z2)) . A very fine spatial mesh 
is therefore required. Since this wave remains localized at the interfaces where it is gen- 
erated, space-time mesh refinement is a good strategy. For this purpose, we adapt a 
steady-state version of the algorithm proposed in ||TH2l- In the fine grid, both the spa- 
tial meshes and the time step are divided by a refinement factor q. Doing so ensures the 
same CFL number in each grid. The coupling between coarse and fine meshes is based 
on spatial and temporal interpolations. 

Even if the refined zone is much smaller than the whole domain, the refinement 
greatly increases the computational cost. The factor q therefore must be chosen ade- 
quately. We choose q that ensures the same discretization of the slow wave, on the refined 
zone, than the fast wave, on the coarse grid. Direct calculations give 

Cps[j0) 

where /o is the central frequency of the signal. 



3.3 Discretization of the interface conditions 

The discretization of the interface conditions requires special care. A straightforward 
stair-step representation of interfaces introduces first-order geometrical errors and yields 
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spurious numerical diffractions. In addition, the jump conditions (|2.10|) are not enforced 
numerically if no special treatment is applied. Lastly, the smoothness requirements to 
solve (I3.2a|) are not satisfied, decreasing the convergence rate of the ADER scheme. 

To remove these drawbacks while maintaining the efficiency of Cartesian grid meth- 
ods, immersed interface methods constitute a good strategy II251I26I . Various formula- 
tions have been proposed in the literature; here, we follow the methodology proposed in 
acoustics / elastic media II27II , viscoelastic media II28II , and poroelastic media IITTH . 
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Figure 3: Irregular point M{xi,yj)ECi-i and its orthogonal projection P onto T. The grid nodes used to compute 
U^j are inside the circle with radius d and centered on P; they are denoted by +. 

To illustrate the basic principle of this method, let us take an irregular point {xj,yj) G 
Qq, where time-marching (I3.2all uses a value at {xj,yj) G Qi (figure |3]). Instead of us- 
ing Uj^j, the discrete operator Ha uses a modified value Uj j. This latter is a r-th order 
extension of the solution from Oq into Qi. 

In a few words, U*j j is build as follows. Let P be the orthogonal projection of {xi,yj) 
on r, and consider the disc T> centered on P with a radius d (figure |3]). Based on the 
interface conditions (|2.30|l at P and on the numerical values U'^' at the grid nodes inside 
V, a matrix Ai is build so that 

U|j = A^(u(i))^. (3.5) 

The derivation of matrix ^A in (|3.5|) is detailed in appendix |Bl Some comments are done 
about the immersed interface method: 

1. A similar algorithm is applied at each irregular point along T and at each propaga- 
tive part of the splitting algorithm (|3.2aD . Since the jump conditions do not vary 
with time, the evaluation of the matrices in (|3.5|) is done during a preprocessing 
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step. Only small matrix-vector products are therefore required at each splitting 
step. After optimization of the computer codes, this additional cost is made negli- 
gible, lower than 1% of the time-marching. 

2. The matrix A4 in (|3.5|) depends on the subcell position of P inside the mesh and on 
the jump conditions at P, involving the local geometry and the curvature of Y at P. 
Consequently, all these insights are incorporated in the modified value (13.511 , and 
hence in the scheme. 

3. The simulations indicate that the number of grid nodes inside the disc V has a 
crucial influence on the stability of the immersed interface method. Here we use a 
constant radius d. Taking r = 2, numerical experiments have shown that d = 3.2Ax 
is a good candidate, while d = 4.5Ax is used when r = 3. 

4. The order r plays an important role on the accuracy of the coupling between the 
immersed interface method and a fc-th order scheme. If r > fc, then a k-th order local 
truncation error is obtained at the irregular points. However, r = k—l suffices to 
keep the global error to the fc-th order ||20| , and hence r=3 is required by the ADER 
4 scheme when S = 0. 

5. A special attention needs to be paid in the case of imperfect hydraulic contact 112.2811 . 
Typical values of JC range around 10^'^ m/s/ Pa. From (|2.30|l and (|A.3|) . it follows 
that numbers close to 10'' coexist with numbers close to 1 in the matrix Lj. The ma- 
trices involved in steps 3 and 4 of the immersed interface method are then badly 
conditioned, which generates numerical instabilities during time-marching. To 
overcome this difficulty, we normalize the physical parameters and the unknowns 
in our codes. This normalization is described in appendix O In practice, it is used 
whatever the interface conditions. 

4 Numerical experiments 
4.1 Physical parameters 

The acoustic medium Qq is water (p/ = 1000 kg/m'^, c = 1500 m/s). The poroelastic 
medium Qi is a water saturated unconsolidated sand, whose material properties are 
summarized in table [T] and are issued from table 3 of Ii37il . In some experiments, an 
inviscid saturating fluid is artificially considered: // = Pa.s, the other parameters being 
unchanged. It is mainly addressed here for a numerical purpose. The cases of open pores 
(|2.26|) . sealed pores l|2.27|l , and imperfect pores (|2.28|) are successively investigated. In the 
latter case, the value of hydraulic permeability /C is given in table [TJ 

Once the spatial mesh sizes Ax = Ay are chosen on the coarse grid, the time step 
follows from the CFL number in Qi: c^.Af/ Ax = 0.95 < 1. The time evolution of the 
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Saturating fluid 


Pf (kg/m^) 


1000 




c (m/ s) 


1500 




f] (Pa.s) 


/ 1.0510-3 


Grain 


Ps (kg/m3) 


2690 




H (Pa) 


1.8610^ 


Matrix 


4> 


0.38 




a 


1.8 




K(m2) 


2.7910-" 




Ao (Pa) 


1.210^ 




m (Pa) 


5.3410'^ 






0.95 


Interface 


/C (m/s/Pa) 


5.10-7 


Phase velocities 


Cp/(/o) (m/s) 


2066.43 




Cps(/o) (m/s) 


124.36 




Cs(/o) (m/s) 


953.05 




c^f (m/s) 


2071.85 




c~ (m/s) 


741.65 




cr (m/s) 


1006.32 




/c (Hz) 


1264.49 



Table 1: Poroelastic medium nj: physical parameters and acoustic properties at /o = 20 Hz. 
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source is a combination of truncated sinusoids 

^ flmSin(/B,„a?oO ifO<f<— , 

m=l /O (4.1) 

otherwise, 

where = 2"'~^, ojq = 27r/o; the coefficients are: fli = 1, fl2 = —21/32, = 63/768, 
fl4= — 1/512, ensuring smoothness of the solution. Two types of sources and boundary 
conditions are considered: 

• no source term fp in (|2.2|l , but an incident plane wave in Qq as initial conditions: 

/ cose \ 



U(x,y,^o) = - 

\ Pf / 

where 6 is the angle between the wavevector and the x-axis, and to adjusts the 
location of the plane wave in Qq. In section l42l the diffracted plane waves are 
computed exactly, and they are enforced numerically on the edges of the compu- 
tational domain. In section l43l periodic computational edges are imposed along 
y-direction, and hence the incident acoustic wave does not need to be enforced; 

• null initial conditions, but a varying source term in (I2.2D 

fp = Six-Xs)5iy-y,)h{t) (4.3) 

that generates cylindrical waves. The size of the domain and the duration of the 
simulations are defined so that no special attention is required to simulate outgoing 
waves, for instance with Perfectly-Matched Layers ||3T]| . 



h{t) = { 



c 

sin0 



xcost 



-ysmt 



(4.2) 



4.2 Test 1: plane wave on a plane interface 

As a first test, we consider a domain [0,400] m^ cut by a plane interface F with slope 60 de- 
grees. An incident plane wave propagates in the fluid, with =0.03 s and 0= —30 degrees 
(14.21) . Consequently, the incident wave crosses the interface normally, leading to a 1-D 
configuration (figure SJ-a); from a numerical point of view, however, the problem is fully 
bidimensional. The advantage of such a 1-D configuration is that each diffracted wave 
has interacted with the interface and is consequently very sensitive to the discretization 
of the interface conditions (12.101) . The central frequency /o = 40 Hz is much smaller than 
the critical Biot frequency: see (|2.5|) and table [H in (|4.1|) is 

We consider an inviscid saturating fluid (?7 = Pa.s): as a consequence, exact solu- 
tions are computed very accurately without Fourier synthesis, and splitting errors of the 



16 
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X (m) X (m) 

Figure 5: test 1. Slices of p at the final instant, for various interface conditions, (a): reflected acoustic wave 
Rp, transmitted fast compressional wave Tpj, transmitted slow compressional wave Tps. (b): zoom around the 
wave reflected in the fluid domain {Rp). 
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scheme are avoided (section [3 .111 . The computations are done on a uniform grid of 400^ 
points, during 150 time steps. Comparisons with the exact values of the pressure p are 
done on the sub-domain [50,350] mx [150,250] m, in order to avoid spurious effects in- 
duced by the edges of the computational domain (figure SJ-b). One observes the reflected 
acoustic wave, the transmitted fast wave and the transmitted slow compressional waves 
(no shear wave is generated in ID). 
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Table 2: test 1. Error measurements and convergence rate in I2 norm, for various interface conditions. Linear 
(r = l), quadratic {r = 2) or cubic ('' = 3) immersed interface method. IVIeasures are done after 3N/8 time 
steps. 



The influence of the interface conditions on the diffracted waves is illustrated in figure 
|5l The exact value of p at the final instant is shown at y = 200 m for the conditions 
(|2.26|) . (|2.27|) and (|2.28|) . The main difference between these three cases is observed in the 
transmitted slow wave. In real experiments, however, only the reflected acoustic wave is 
measured: a zoom on this wave is given in figure |5l-b. The hydraulic permeability /C = 
5.10-^^ leads to an intermediate regime between open and sealed pores. If /C > 10-^, the 
results (not shown here) cannot be distinguished from those obtained with open pores. 
In the same way, results obtained with /C < 10-^ cannot be distinguished from the sealed 
pores. 

How accurate is the discretization of the interface conditions is assessed through com- 
parisons with the analytical solutions (figure[6|. In the left row, such comparisons are pro- 
posed in two cases; r = means that no numerical treatment is done along the interface, 
and the numerical solution does not converge towards the exact one. On the contrary, 
excellent agreement is observed when r = 2 is used (section l331l . 
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Figure 6: test 1. Left row: comparisons between exact values and numerical values of p. Right row: errors 
measured in I2 norm versus the number of grid nodes, with various order r of the immersed interface method. 
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Error measurements on successive refined grids are given in table |2l Convergence 
rates are drawn on the right row of figure [6l Various values of the order r of the im- 
mersed interface method are investigated. As stated in section l33l fourth-order accu- 
racy is maintained if third-order extrapolations (r = 3) are used in the immersed interface 
method. From now on, all the simulations are performed with r = 3. 



(a) 



(b) 




100 150 200 250 300 350 



analytical 

o numerical q=6 















X (m) 



(c) 



(d) 






analytical 








'J numerical q=6 








numerical q=2 






5 


<- numerical q=1 






Pressure (Pa) 






-2 







X (m) 



X (m) 



Figure 7: test 1, with a viscous saturating fluid (;/ = 1.0510^^ Pa.s) and open pore conditions. Snapshot (a) 
and slice (b) of the exact value of p at the final instant; zooms on the reflected acoustic wave (c) and on the 
transmitted slow compressional wave (d), for various grid sizes. Vertical line represents the location of the fluid 
/ porous interface. 

In figurelZl we take into account the viscosity of the the saturating fluid (rj = 1.0510^^ 
Pa.s), the other parameters being unchanged (open pore conditions are considered). In 
this case, the slow compressional wave is a static mode that remains where it is generated, 
close to the interface in the porous medium, and hence a much finer grid is required 
(section |32|) . The exact solution is computed by Fourier synthesis on 256 modes, with a 
frequency step A/ = 1.6 Hz. 

In comparison with figures |4l-(b) and|6l-(a), the transmitted slow compressional wave 
cannot be seen in figure [7l-(b). Zooms around the reflected acoustic wave (c) and around 
the transmitted slow compressional wave (d) are proposed for various grids. If the coarse 
grid of the inviscid case is used (q = 1), then large errors are observed compared with the 
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exact solution. It is a consequence of the too crude discretization of small-scale slow wave 
(d). On the contrary, if a sufficiently fine grid is used {q = 6), an excellent agreement is 
obtain between numerical and analytical solutions. 

4.3 Test 2: plane wave on a circular interface 

A circular interface T of radius 100 m is centered at (0,0) in the domain [— 600,600] m^ 
(figure [8]). This configuration is relevant to examine the discretization of an interface with 
constant non-zero curvature. The cylinder is filled with the porous medium Qi, while 
the acoustic medium Qq lies outside. The computational domain is discretized on 800^ 
points, leading to Ax = Ay = 1.5 m. A locally refined mesh [— 110,110] m^ is used around 
the interface in order to compute accurately the different wave conversions. As in test 1, 
the source is an acoustic plane wave initially in Qq, with 6 = degree (|4.2|) . The central 
frequency in (14.11) is /o = 20 Hz. The initial conditions are illustrated in figure [51 



(a) (b) 




Figure 8: test 2. Initial values of p: snapshot (a) and slice along the x-axis, at i/ = m (b). The dotted square 
in (a) represents the frontiers of the refined mesh. The vertical lines in (b) denote the location of interfaces. 

The first experiments concerns the inviscid case (fj = 0), where the poroelastic waves 
are non dispersive. Based on the criterion (|3.4|l and on table[TJ the refinement factor in the 
local grid is set to = 3 « / c~ . The numerical values of p after 430 time steps (t ~ 0.22 
s) are displayed in the left row of figure |9l for the three types of interface conditions. The 
pressure recorded during 900 time steps at RO {x,- = — 120,i/r = 0) in Oq is shown in the 
right row for t > 0.2 s: the incident wave and the first refracted wave are not represented, 
in order to focus on the successive reflected /transmitted waves which strongly depend 
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Figure 9: test 2, inviscid saturating fluid (f/ = Pa.s). Snapshots of p after 430 time steps (left row), time 
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on the interface conditions. This is particularly true for t > 0.4 s, where shape and ampli- 
tude of the recorded pressure completely differ depending on the hydraulic contact. 

In the inviscid case, the analytical solutions can be computed very accurately, by a 
decomposition of Fourier modes on a basis of circular functions. In practice, reference 
values are obtained by using Nf =32768 Fourier modes (with a frequency step A/=0.0063 
Hz) and a truncated basis of 70 Bessel functions. The agreement between numerical and 
exact values is excellent in the three cases (figure |9l right row). 



(a) (b) 




■Wa -400 ■20C 200 400 500 



t(s) 

Figure 10: test 2, viscous saturating fluid, (a): snapshots of p after 430 iterations for imperfect pore condition, 
(b): time history of p at the receiver RO in the acoustic medium, for various refinement factors q. 

Similar experiments are also performed with a viscous saturating fluid. The numer- 
ical values of p after 430 time steps are displayed in figure [TOl-a, with imperfect pores. 
From the criterion (|3.4|) and the phase velocities given in table [H one obtains cj =16, which 
is very costly. Nevertheless, the pressure recorded at RO reveals that q = 9 suffices to get 
reference solutions (figure ITOl-b). 

Compared with the inviscid vase, the viscosity greatly modifies the signal recorded at 
receiver RO (figure ITTI-a). Figure [TTl-b shows the reflected waves obtained with the three 
pore conditions. The differences between these signals are smaller than in the inviscid 
case (right row of figure |9]l. 

4.4 Test 3: a sinusoidal interface 

As a third and last test, we consider a sinusoidal interface separating the acoustic medium 
Oq (top) and the poroelastic medium Qi (bottom). This configuration may model the sea- 
floor; numerically, it is relevant to test the algorithms with a non-constant curvature of 
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the interface. The domain is [—1500, 1500] m^, and the interface is given by the relation 
y = 40sin( . The viscosity of the saturating fluid is taken into account (?/ 7^ 0). A 
source term fp is put at the point ( x, = 0, i/s = 20 ) in Oq (I4.3I) . The mesh size is A x = A 1/ = 2 
m, except in the vicinity of the interface, where a refinement factor oiq=7 is applied. The 
refined grid contains about 3.510^ nodes and 75544 irregular nodes where the immersed 
interface method is applied (section |33)| . 

The pressure field after 800 iterations {t ~ 0.73 s) is displayed in figure [121 From the 
criterion (13.411 and the phase velocities given in table [TJ one obtains q = 16, which is very 
costly. Nevertheless, the pressure recorded close to the interface, at Rl (x=750,i/= —45.2) 
in Oi, reveals that grid convergence is satisfying when = 7 (figure [T3]-a). 

The time history of the pressure recorded at RO (x = 500,1/ = 200) in Qq is displayed 
in figure [l3]-b for the three pore conditions. As in test 2 (figure [TTI-b), the first reflected 
waves almost do not depend on the hydraulic contact. Consequently, we focus on the 
subsequent waves (^ > 0.435 s), where differences are clearly observed depending on the 
interface conditions. 



5 Conclusion 

We have developed a robust and highly accurate numerical model to simulate wave 
propagation in fluid / poroelastic media. Our model can incorporate various models of 
interface conditions, in particular the case of imperfect hydraulic contact: to our knowl- 
edge, it is the first time that such simulations are proposed. Arbitrary-shaped interfaces 
can be handled and accuracy is ensured by a subcell resolution on a Cartesian grid. 

Numerical experiments have shown that each part of the algorithm is required to 
get efficiently reliable results: fourth-order scheme with time splitting, mesh refinement. 
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Figure 12: test 3. Snapshot of p after 800 iterations, taking imperfect pore conditions. Horizontal black dotted 
lines denote the frontiers of the refined grid. 



immersed interface method. The effect of interface conditions on the diffracted waves has 
been illustrated. When the viscous effects are noticeable, the accurate computation of the 
slow compressional wave in the poroelastic medium is crucial for the overall accuracy. 
This diffusive wave propagates along the interface and plays a major role in the balance 



25 



(a) 



(b) 





A 




numerical q^3 








numerical q=5 








numerical q=7 





















t(s) 




Figure 13: 
of p at the 



test 3. (a): time history of p at the receiver Rl in Qj for different refinement factors; time history 
receiver RO in CIq, for t> 0.435 s (b). 



equations. 



This numerical model enables to investigate many physically-relevant configurations. 
For instance, comparisons between real experiments and simulations could be used to 
characterize the hydraulic permeability fC in (|2.10|) . Simulation of multiple scattering in 
random or periodic media is another fruitful application. The objective is to estimate 
numerically the properties of the homogenized effective medium p30^. Since a Cartesian 
grid is used and no meshing of the interfaces is required, our approach is very well suited 
to the modeling of numerous scatterers (typically, a few hundreds UTOH '). 



Our numerical model is valid only in the low-frequency range. For frequencies greater 
than fc in (|2.5|) , the second equation of (|2.3|) must be modified to account for the viscous 
layer dissipation. Fractional derivatives of order 1/2 are then introduced ||29|. We are 
currently investigating this topic [Sj. Similar extensions are also required concerning the 
hydraulic permeability in (12.101) . Indeed, it is known that fC depends not only on the 
geometrical properties of the medium, but also on the frequency Il36ll . 
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A Matrices of interface conditions 

The matrices Cq, C'j' and introduced in (|2.30|) are detailed. In the case of open pore 
conditions (|2.26ll , it follows from (I2.1I) 
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In the case of sealed pore conditions (I2.27D , one gets 
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In the case of imperfect pore conditions (|2.28|) , one gets 
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(A.3) 



B Four steps to build U^^ (|33) 

Step 1: high-order interface conditions. The conditions (I2.30D are differentiated in terms 
of t and T. Time derivatives are replaced by spatial derivatives thanks to the splitted 
evolution equations (equation (l3.2aD ). Consider for instance the zero-th order boundary 
conditions L^U^ = 0; time derivative yields 

^l!u;^l»^u^l;(-a4u;-b,^u?)^o. (b.d 
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Derivatives in terms of T are replaced by spatial derivatives thanks to the chain-rule. For 
instance, = yields 

^uu^(^L»)u;+Lo(4u;+y^u;).o. (b.2) 

From (|2.30|) , jB.l^ and (|B.2|) , we build a matrix L\ such that L| U| =0, which provides first- 
order boundary conditions. A similar procedure is applied to the zeroth-th order jump 
conditions. By iterating this process r times, r-th order interface conditions are obtained 

CiU'i = QU^, KU[ = 0. (B.3) 

In our codes, the matrices and L'^ are computed automatically by developing com- 
puter algebra tools. 

Step 2: high-order compatibility conditions. Some components of the successive spatial 
derivatives of U are not independent. In the fluid medium Qq, the vorticity of acoustic 
velocity is null V A v = 0, which yields 

dx dy 

In the porous medium Oi, the symmetry of a yields the Beltrami-Michell equation 

dxdy dx^ dx'^ dx^ dy^ dy^ dy^' 

4(Ao + h)' ' 4(Ao + ?/)' ' 2(Ao + /y)- 

Equations (IB.4II and dB.SIl can be differentiated in terms of x and y as many times as re- 
quired, assuming a sufficiently smooth solution. The equations so-obtained are satisfied 
at any point, in particular on both sides of I . They can be used to reduce the number of 
components of U-'. Reduced vectors V, are defined, such that 

U^G^V^, (B.6) 

where G[ are x (riy — ni,) matrices. In Qq, n;, = r (r+l)/2 if r > 1, else; an algorithm to 
compute Gq is given in appendix A of II27II . In Oi, n;, = r(r — 1)/2 if r >2, nf, = else; an 
algorithm to compute G[ is given in appendix B of IfTTH . The compatibility conditions are 
crucial to ensure the stability of the immersed interface method. 

Step 3: high-order boundary values. On the poroelastic Oi side, the r-th order bound- 
ary conditions in (|B.3|) and the high-order Beltrami equations (|B.6|) yield the underdeter- 
mined linear system 

L[G[V[ = 0. (B.7) 
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It leads to 

Y[ = K[W[, (B.8) 

where is the matrix built from the kernel of l-[G[. The solution W[ is the minimum 
set of independent components of the trace of U and its spatial derivatives up to the r-th 
order, on the domain Oi . Injecting (|B.8|l into the r-th order jump conditions (|B.3|I gives 

S[W[ = S'oyl (B.9) 

where S\ = C[ G[ K.[ and Sq = CqGq. The underdetermined system (|B.9|l is solved 

Wl=((SD"'S^|Rs;)(^^°^, (B.IO) 

where (S![)^^ is the least-squares pseudo-inverse of Rgj is the matrix containing the 
kernel of and A'' is a set of Lagrange multipliers. To build (SJ)^^ and Rsj, a singular 
value decomposition of S\ is performed. 

Step 4: construction of modified values. The coefficients of 2-D Taylor expansions 
around P (figure [S]! are put in the matrix IV- ■: 

where a=0,...,r and /3 = 0, . . ., a; I is the 3 x 3 or 8 x 8 identity matrix, depending on whether 
{xi,yj) belongs to the fluid domain Qq or the poroelastic domain Oi. By definition, the 
modified value at (x/,!//) is 

To determine the trace Uq in (IB.12II , consider the disc T> centered at P with radius d (figure 
|3]). At the grid nodes of PnOo, r-th order Taylor expansions at P, (|B.6|) and (|B.8|) give 

t(1) 



Uj% = n^U^. (B.12) 




(B.13) 



where Taylor rests have been omitted for the sake of clarity. At the grid nodes of "DnOi, 
r-th order Taylor expansions at P, (|B.6|l , (|B.8|l and (|B.10|) . give 

111' ^ ^ ^ (B.14) 

= n;:.GiKi((SO-'|Rs;^ 





I; 


|Rs;) 
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The equations (|B.13|) and (|B.14|) are written using an adequate matrix M 

The least-squares inverse of the matrix M is denoted by M^^. Since the Lagrange multi- 
pliers A*" are not involved in (IB.12I) , M^^ is restricted to M ^, so that 

V(j = M"^fuW) . (B.16) 



The modified value follows from (|B.6|) , (|B.8|) , (|B.12|) and (|B.16|) , recovering (|3.5|) : 



V 

(B.17) 

C Normalization of the variables 

As mentioned in section |331 normalized parameters and unknowns are used in our com- 
puter codes. These quantities are denoted by tildes in the following. Given a real M , we 
define the normalized time 

t=Mt, (C.l) 

the normalized variables 

v = Arv, \s=M\s, ^=M^N, ^=^r P = Jf' (C-2) 
and the normalized physical parameters 

~ _Pf ^ _Ps _ _Piv 

A~=il u = JL rh = — (C.3) 

The value of the normalization parameter is set to A/" = 1000 in numerical experiments. 
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